####################################################
### Estimate Demand Model for Covered California ###
####################################################
	
#ssubmit --partition=sscc --cores=10 --mem=200g "julia18 calculate.model.fit.holdout.forward2.jl"		
	
cluster_run = true;	
	
if cluster_run
	codedir = "/project/inertia_project/Code";
	datadir = "/project/inertia_project/Data";
else
	codedir = "C:/Users/eas24f/OneDrive - Florida State University/Research/Inertia RESTAT/Code"; #"V:/inertia_project/Code";
	datadir = "C:/Users/eas24f/OneDrive - Florida State University/Research/Inertia RESTAT/Data"; #"V:/inertia_project/Data";
end
		
cd(datadir)
using DataFrames
using StatsBase
using Optim
using Calculus
using Distributions
using Random
using DelimitedFiles
using CSV
using LinearAlgebra		
		
##### Specification flags
	
	add_premium_interactions = true;
	add_previous_demographics_interactions = true;
	add_previous_firm_interactions = true;
	add_av_interactions = true;
	add_intercepts = true;
	add_rating_areas = true;
	add_insurer_market_fe = true;
	add_silver = false;
	add_complex_interactions = true;
	
	control_function = true;
	use_eta = false;
	random_parameters = true;
	sequence_flag = true;
	inattention = false;
	nested = false;
	ns = 50;
	
	year_to_run = 2018; # supply model uses data from 2014-2018
	
	churner_run = false;
	lapser_run = false;
	age_income_run = false;
	
	rich_flag = false;
	poor_flag = false;
	old_flag = false;
	young_flag = false;
	lapse_flag = false;
	no_lapse_flag = false;
	churn_flag = false;
	no_churn_flag = false;
	eq_robust_flag = false;
	
	forward_flag = false;
	
##### Load Data 
	
	data = CSV.read("fit_julia_data.csv",DataFrame,header=true);
	households = CSV.read("fit_households.csv",DataFrame,header=true); # household characteristics 
	plans = CSV.read("ca_plan_characteristics_AUG032019_small.csv",DataFrame,header=true); # plan characteristics (does not include same plans in different rating areas and years)
	param_estimates = CSV.read("demand_parameter_estimates.csv",DataFrame,header=true); # parameters from demand estimation
	Random.seed!(6);
	
	choices = convert(Array,data[!,:choice]);
	previous_choices = convert(Array,data[!,:previous_choice]);
	weights = convert(Array,data[!,:weight]);
	uninsured_plans = convert(Array,data[!,:uninsured_plan]);
	household_numbers = convert(Array,data[!,:household_number]);			
	
#### Read in Parameters 


	if add_av_interactions & add_insurer_market_fe & sequence_flag & random_parameters & control_function
		param_estimates = dropmissing(param_estimates,:AV_Seq_CF_Random_IMFE);
		min_beta = param_estimates[!,:AV_Seq_CF_Random_IMFE];
	elseif add_av_interactions & add_insurer_market_fe & random_parameters & control_function
		param_estimates = dropmissing(param_estimates,:AV_CF_Random_IMFE);
		min_beta = param_estimates[!,:AV_CF_Random_IMFE];
	elseif add_av_interactions & sequence_flag & random_parameters & control_function
		param_estimates = dropmissing(param_estimates,:AV_Seq_CF_Random);
		min_beta = param_estimates[!,:AV_Seq_CF_Random];
	#elseif inattention
	#	param_estimates = dropmissing(param_estimates,:Inattention);
	#	min_beta = param_estimates[!,:Inattention];
	#elseif sequence_flag & random_parameters & control_function & add_insurer_market_fe
	elseif sequence_flag & random_parameters & control_function & add_insurer_market_fe
		param_estimates = dropmissing(param_estimates,:Seq_CF_Random_IMFE);
		min_beta = param_estimates[!,:Seq_CF_Random_IMFE];
	elseif random_parameters & control_function & add_insurer_market_fe
		param_estimates = dropmissing(param_estimates,:CF_Random_IMFE);
		min_beta = param_estimates[!,:CF_Random_IMFE];
	elseif random_parameters & control_function & sequence_flag & add_silver
		param_estimates = dropmissing(param_estimates,:Silver_Seq_CF_Random);
		min_beta = param_estimates[!,:Silver_Seq_CF_Random];
	elseif random_parameters & control_function & sequence_flag
		param_estimates = dropmissing(param_estimates,:Seq_CF_Random);
		min_beta = param_estimates[!,:Seq_CF_Random];
	elseif random_parameters & control_function & add_silver
		param_estimates = dropmissing(param_estimates,:Silver_CF_Random);
		min_beta = param_estimates[!,:Silver_CF_Random];
	elseif random_parameters & control_function
		param_estimates = dropmissing(param_estimates,:CF_Random);
		min_beta = param_estimates[!,:CF_Random];
	elseif random_parameters & sequence_flag
		param_estimates = dropmissing(param_estimates,:Seq_Random);
		min_beta = param_estimates[!,:Seq_Random];
	#elseif random_parameters
	#	param_estimates = dropmissing(param_estimates,:Random);
	#	min_beta = param_estimates[!,:Random];
	#elseif nested
	#	param_estimates = dropmissing(param_estimates,:Nested_Logit);
	#	min_beta = param_estimates[!,:Nested_Logit];
	else
		param_estimates = dropmissing(param_estimates,:Logit);
		min_beta = param_estimates[!,:Logit];
	end

	if age_income_run
		param_estimates = dropmissing(param_estimates,:RichOld);
	elseif lapser_run
		param_estimates = dropmissing(param_estimates,:Lapser);
	elseif churner_run
		param_estimates = dropmissing(param_estimates,:Churn);
	end
		
	
##### Specify covariates
	
	all_covariates = convert(Array,param_estimates[!,:Variable]);
	if nested
		pop!(all_covariates);
	end
	if inattention
		premium_shock_index = findall(all_covariates .== "premium_shock")[1];
		search_covariates = all_covariates[premium_shock_index:length(all_covariates)];
		all_covariates = all_covariates[1:(premium_shock_index-1)];
	end
	
	if add_complex_interactions & add_silver
		product_covariates = ["Premium","AV","Silver","HMO","Anthem","Blue_Shield","Kaiser","Health_Net"];	
	elseif add_complex_interactions
		product_covariates = ["Premium","AV","HMO","Anthem","Blue_Shield","Kaiser","Health_Net"];	
	elseif add_silver
		product_covariates = ["Premium","AV","Silver","HMO","Anthem","Blue_Shield","Kaiser","Health_Net","Anthem_HMO"];	
	else
		product_covariates = ["Premium","AV","HMO","Anthem","Blue_Shield","Kaiser","Health_Net","Anthem_HMO"];	
	end
	
	if random_parameters
		if add_complex_interactions & add_av_interactions
			random_variables = ["intercept","Anthem","Blue_Shield","Kaiser","Health_Net"];
			random_covariance_variables = ["Anthem_random","BS_Anthem_cov","Kaiser_Anthem_cov","HN_Anthem_cov","AV_Anthem_cov",
											"Blue_Shield_random","Kaiser_BS_cov","HN_BS_cov","AV_BS_cov","Kaiser_random","HN_Kaiser_cov",
											"AV_Kaiser_cov","Health_Net_random","AV_HN_cov","AV_random"];
			random_parameter_product_indices = cat(dims=1,findall(all_covariates .== "intercept")[1],findall(in(random_variables),product_covariates));
		elseif add_complex_interactions 
			random_variables = ["intercept","Anthem","Blue_Shield","Kaiser","Health_Net","AV"];
			random_covariance_variables = ["Anthem_random","BS_Anthem_cov","Kaiser_Anthem_cov","HN_Anthem_cov","AV_Anthem_cov",
											"Blue_Shield_random","Kaiser_BS_cov","HN_BS_cov","AV_BS_cov","Kaiser_random","HN_Kaiser_cov",
											"AV_Kaiser_cov","Health_Net_random","AV_HN_cov","AV_random"];
			random_parameter_product_indices = cat(dims=1,findall(all_covariates .== "intercept")[1],findall(in(random_variables),product_covariates));
		elseif inattention & add_silver
			random_variables = ["AV","Silver","HMO","Anthem","Blue_Shield","Kaiser","Health_Net"];
			random_parameter_product_indices = findall(in(random_variables),product_covariates);
		elseif inattention
			random_variables = ["AV","HMO","Anthem","Blue_Shield","Kaiser","Health_Net"];
			random_parameter_product_indices = findall(in(random_variables),product_covariates);
		elseif add_silver
			random_variables = ["intercept","AV","Silver","HMO","Anthem","Blue_Shield","Kaiser","Health_Net"];
			random_parameter_product_indices = cat(dims=1,findall(all_covariates .== "intercept")[1],findall(in(random_variables),product_covariates));
		else
			random_variables = ["intercept","AV","HMO","Anthem","Blue_Shield","Kaiser","Health_Net"];
			random_parameter_product_indices = cat(dims=1,findall(all_covariates .== "intercept")[1],findall(in(random_variables),product_covariates));
		end
		random_product_covariates = string.(random_variables,"_random");
		
		random_parameter_indices = findall(in(random_product_covariates),all_covariates);
		if length(random_parameter_indices) == 1
			random_parameter_indices = random_parameter_indices[1];
		end
	end

##### Functions 

	function simulate_nested_probabilities(betap,Xins,household_numbers_ins,uninsured_plans,indices) 
		lambda = last(betap);
		betap = betap[1:(length(betap)-1)];

		expdf = DataFrame(household_number = household_numbers_ins,expvarall = exp.((Xins * betap)/lambda));
		expsum = combine(groupby(expdf,:household_number), :expvarall => sum)[!,:expvarall_sum];
		
		probs = zeros(length(uninsured_plans));
		probs[uninsured_plans .== 0] = expdf[!,:expvarall] .* (expsum[household_numbers_ins]).^(lambda-1) ./ (1 .+ (expsum[household_numbers_ins]).^(lambda));
		probs[uninsured_plans .== 1] = 1  ./ (1 .+ (expsum).^(lambda));
		
		return(probs[indices])
	end		
	
	function simulate_logit_probabilities(betap,Xins,household_numbers_ins,uninsured_plans,indices) 
		expdf = DataFrame(household_number = household_numbers_ins,expvarall = exp.(Xins * betap));
		expsum = combine(groupby(expdf,:household_number), :expvarall => sum)[!,:expvarall_sum];
		
		probs = zeros(length(uninsured_plans));
		probs[uninsured_plans .== 0] = expdf[!,:expvarall] ./ (1 .+ expsum[household_numbers_ins]);
		probs[uninsured_plans .== 1] = 1  ./ (1 .+ expsum);
		
		return(probs)
	end	

	function get_choice_probabilities(min_beta,X,V,choices,household_numbers,uninsured_plans)
		
		Xupdated = copy(X);
		if nested
			probabilities = simulate_nested_probabilities(min_beta,X[uninsured_plans .== 0,:],household_numbers[uninsured_plans .== 0],uninsured_plans,collect(1:1:length(choices)));
		elseif inattention
			probabilities = zeros(length(choices));
			for n in 1:ns
				Xupdated[:,random_parameter_indices] = X[:,random_parameter_indices] .* V[:,:,n]; 
				probabilities += simulate_logit_probabilities(min_beta[1:(premium_shock_index-1)],Xupdated[uninsured_plans .== 0,:],household_numbers[uninsured_plans .== 0],uninsured_plans,collect(1:1:length(probabilities))) ./ ns;
			end
		elseif random_parameters
			probabilities = zeros(length(choices));
			for n in 1:ns
				if control_function & use_eta
					Xupdated[:,findall(all_covariates .== "eta")[1]] = eta[:,n];
				end
				Xupdated[:,random_parameter_indices] = X[:,random_parameter_indices] .* V[:,:,n]; 
				probabilities += simulate_logit_probabilities(min_beta,Xupdated[uninsured_plans .== 0,:],household_numbers[uninsured_plans .== 0],uninsured_plans,collect(1:1:length(probabilities))) ./ ns;
			end
		else
			probabilities = simulate_logit_probabilities(min_beta,X[uninsured_plans .== 0,:],household_numbers[uninsured_plans .== 0],uninsured_plans,collect(1:1:length(choices)));
		end
		
		return(probabilities)
	end
	
	function get_dominated_enrollment(simulated_choice_weights,X,data,households,household_numbers)
		gold_dominated = sum(simulated_choice_weights[intersect(findall(data[!,:av] .== 0.8),findall(households[household_numbers,:FPL] .<= 2))]);
		platinum_dominated = sum(simulated_choice_weights[intersect(findall(data[!,:av] .== 0.9),findall(households[household_numbers,:FPL] .<= 1.5))]);
		return(cat(dims=1,0,gold_dominated + platinum_dominated,-gold_dominated,-platinum_dominated))
	end
	
	# Form X covariate matrix
	function form_X_matrix(all_covariates,data,households,choices,previous_choices,weights,uninsured_plans,household_numbers,V)	
	
		prem_ind = findall(all_covariates .== "Premium")[1];
		if !inattention
			prev_ind = findall(all_covariates .== "previous_choice")[1];
			#silver_ind = findall(all_covariates .== "Silver")[1];
		end
		av_ind = findall(all_covariates .== "AV")[1];
		
		X = zeros(Float64,length(choices),length(all_covariates));
		X[:,prem_ind] = convert(Array,data[!,:premium]) ./ households[household_numbers,:enrollees];
		if !inattention
			X[:,prev_ind] = convert(Array,data[!,:previous_choice]);
		end
		X[:,findall(all_covariates .== "AV")[1]] = convert(Array,data[!,:av]);
		X[:,findall(all_covariates .== "intercept")] = 1 .- uninsured_plans;		
		
		# Add product characteristics
		plans[!,:PPO] = plans[!,:HMO];
		for plan in plans[!,:Plan_Name]
			plan_indices = findall(data[!,:plan_name] .== plan);
			for covariate in setdiff(product_covariates,["Premium","AV","Anthem_HMO"])
				X[plan_indices,findall(all_covariates .== covariate)[1]] .= 
					plans[findall(plans[!,:Plan_Name] .== plan)[1],findall(names(plans) .== covariate)[1]];
			end
		end
		
		if !add_complex_interactions
			X[:,findall(all_covariates .== "Anthem_HMO")[1]] = X[:,findall(all_covariates .== "Anthem")[1]] .* X[:,findall(all_covariates .== "HMO")[1]];
		end
		
		# Add premium interactions
		if add_premium_interactions
			if (rich_flag | poor_flag | old_flag | young_flag)
				X[:,findall(all_covariates .== "Premium_male")[1]] = X[:,prem_ind] .* (households[household_numbers,:perc_male]);
				X[:,findall(all_covariates .== "Premium_family")[1]] = X[:,prem_ind] .* (households[household_numbers,:household_size] .> 1);
				X[:,findall(all_covariates .== "Premium_Asian")[1]] = X[:,prem_ind] .* (households[household_numbers,:perc_asian]);
				X[:,findall(all_covariates .== "Premium_Black")[1]] = X[:,prem_ind] .* (households[household_numbers,:perc_black]);
				X[:,findall(all_covariates .== "Premium_Hispanic")[1]] = X[:,prem_ind] .* (households[household_numbers,:perc_hispanic]);
				X[:,findall(all_covariates .== "Premium_Other")[1]] = X[:,prem_ind] .* (households[household_numbers,:perc_other]);
			else
				X[:,findall(all_covariates .== "Premium_250to400")[1]] = X[:,prem_ind] .* ((households[household_numbers,:FPL] .> 2.50) .& (households[household_numbers,:FPL] .<= 4)) ; 
				X[:,findall(all_covariates .== "Premium_gt400")[1]] = X[:,prem_ind] .* (households[household_numbers,:FPL] .> 4) ; 
				if add_complex_interactions
					X[:,findall(all_covariates .== "Premium_0to34")[1]] = X[:,prem_ind] .* (households[household_numbers,:perc_0to17] .+ 
						households[household_numbers,:perc_18to25] .+ households[household_numbers,:perc_26to34]);
				else
					X[:,findall(all_covariates .== "Premium_0to17")[1]] = X[:,prem_ind] .* (households[household_numbers,:perc_0to17]);
					X[:,findall(all_covariates .== "Premium_18to34")[1]] = X[:,prem_ind] .* (households[household_numbers,:perc_18to25] .+ households[household_numbers,:perc_26to34]);
				end
				X[:,findall(all_covariates .== "Premium_35to54")[1]] =  X[:,prem_ind] .* (households[household_numbers,:perc_35to44] .+ households[household_numbers,:perc_45to54]);
				X[:,findall(all_covariates .== "Premium_male")[1]] = X[:,prem_ind] .* (households[household_numbers,:perc_male]);
				X[:,findall(all_covariates .== "Premium_family")[1]] = X[:,prem_ind] .* (households[household_numbers,:household_size] .> 1);
				X[:,findall(all_covariates .== "Premium_Asian")[1]] = X[:,prem_ind] .* (households[household_numbers,:perc_asian]);
				X[:,findall(all_covariates .== "Premium_Black")[1]] = X[:,prem_ind] .* (households[household_numbers,:perc_black]);
				X[:,findall(all_covariates .== "Premium_Hispanic")[1]] = X[:,prem_ind] .* (households[household_numbers,:perc_hispanic]);
				X[:,findall(all_covariates .== "Premium_Other")[1]] = X[:,prem_ind] .* (households[household_numbers,:perc_other]);
			
				if add_complex_interactions
					X[:,findall(all_covariates .== "Premium_250to400_0to34")[1]] = X[:,prem_ind] .* ((households[household_numbers,:FPL] .> 2.50) .& (households[household_numbers,:FPL] .<= 4)) .*
						(households[household_numbers,:perc_0to17] .+ households[household_numbers,:perc_18to25] .+ households[household_numbers,:perc_26to34]); 
					X[:,findall(all_covariates .== "Premium_gt400_0to34")[1]] = X[:,prem_ind] .* (households[household_numbers,:FPL] .> 4) .*
						(households[household_numbers,:perc_0to17] .+ households[household_numbers,:perc_18to25] .+ households[household_numbers,:perc_26to34]);  
					X[:,findall(all_covariates .== "Premium_250to400_35to54")[1]] = X[:,prem_ind] .* ((households[household_numbers,:FPL] .> 2.50) .& (households[household_numbers,:FPL] .<= 4)) .*
						(households[household_numbers,:perc_35to44] .+ households[household_numbers,:perc_45to54]); 
					X[:,findall(all_covariates .== "Premium_gt400_35to54")[1]] = X[:,prem_ind] .* (households[household_numbers,:FPL] .> 4) .*
						(households[household_numbers,:perc_35to44] .+ households[household_numbers,:perc_45to54]);  
				end
			end
			
		end
	
		# AV interactions
		if add_av_interactions
			X[:,findall(all_covariates .== "AV_250to400")[1]] = X[:,av_ind] .* ((households[household_numbers,:FPL] .> 2.50) .& (households[household_numbers,:FPL] .<= 4)) ; 
			X[:,findall(all_covariates .== "AV_gt400")[1]] = X[:,av_ind] .* (households[household_numbers,:FPL] .> 4) ; 
			X[:,findall(all_covariates .== "AV_0to17")[1]] = X[:,av_ind] .* (households[household_numbers,:perc_0to17]);
			X[:,findall(all_covariates .== "AV_18to34")[1]] = X[:,av_ind] .* (households[household_numbers,:perc_18to25] .+ households[household_numbers,:perc_26to34]);
			X[:,findall(all_covariates .== "AV_35to54")[1]] =  X[:,av_ind] .* (households[household_numbers,:perc_35to44] .+ households[household_numbers,:perc_45to54]);
			X[:,findall(all_covariates .== "AV_male")[1]] = X[:,av_ind] .* (households[household_numbers,:perc_male]);
			X[:,findall(all_covariates .== "AV_family")[1]] = X[:,av_ind] .* (households[household_numbers,:household_size] .> 1);
			X[:,findall(all_covariates .== "AV_Asian")[1]] = X[:,av_ind] .* (households[household_numbers,:perc_asian]);
			X[:,findall(all_covariates .== "AV_Black")[1]] = X[:,av_ind] .* (households[household_numbers,:perc_black]);
			X[:,findall(all_covariates .== "AV_Hispanic")[1]] = X[:,av_ind] .* (households[household_numbers,:perc_hispanic]);
			X[:,findall(all_covariates .== "AV_Other")[1]] = X[:,av_ind] .* (households[household_numbers,:perc_other]);
			
			if add_complex_interactions
				X[:,findall(all_covariates .== "AV_250to400_0to34")[1]] = X[:,av_ind] .* ((households[household_numbers,:FPL] .> 2.50) .& (households[household_numbers,:FPL] .<= 4)) .*
					(households[household_numbers,:perc_0to17] .+ households[household_numbers,:perc_18to25] .+ households[household_numbers,:perc_26to34]); 
				X[:,findall(all_covariates .== "AV_gt400_0to34")[1]] = X[:,av_ind] .* (households[household_numbers,:FPL] .> 4) .*
					(households[household_numbers,:perc_0to17] .+ households[household_numbers,:perc_18to25] .+ households[household_numbers,:perc_26to34]);  
				X[:,findall(all_covariates .== "AV_250to400_35to54")[1]] = X[:,av_ind] .* ((households[household_numbers,:FPL] .> 2.50) .& (households[household_numbers,:FPL] .<= 4)) .*
					(households[household_numbers,:perc_35to44] .+ households[household_numbers,:perc_45to54]); 
				X[:,findall(all_covariates .== "AV_gt400_35to54")[1]] = X[:,av_ind] .* (households[household_numbers,:FPL] .> 4) .*
					(households[household_numbers,:perc_35to44] .+ households[household_numbers,:perc_45to54]);  
			end
		end
	
	
		# Add inertia interactions
		if add_previous_demographics_interactions & !inattention
			if (rich_flag | poor_flag | old_flag | young_flag)
				X[:,findall(all_covariates .== "previous_choice_male")[1]] = X[:,findall(all_covariates .== "previous_choice")[1]] .* (households[household_numbers,:perc_male]);
				X[:,findall(all_covariates .== "previous_choice_family")[1]] = X[:,findall(all_covariates .== "previous_choice")[1]] .* (households[household_numbers,:household_size] .> 1);
				X[:,findall(all_covariates .== "previous_choice_Asian")[1]] = X[:,findall(all_covariates .== "previous_choice")[1]] .* (households[household_numbers,:perc_asian]);
				X[:,findall(all_covariates .== "previous_choice_Black")[1]] = X[:,findall(all_covariates .== "previous_choice")[1]] .* (households[household_numbers,:perc_black]);
				X[:,findall(all_covariates .== "previous_choice_Hispanic")[1]] = X[:,findall(all_covariates .== "previous_choice")[1]] .* (households[household_numbers,:perc_hispanic]);
				X[:,findall(all_covariates .== "previous_choice_Other")[1]] = X[:,findall(all_covariates .== "previous_choice")[1]] .* (households[household_numbers,:perc_other]);
			else
				X[:,findall(all_covariates .== "previous_choice_250to400")[1]] = X[:,findall(all_covariates .== "previous_choice")[1]] .* ((households[household_numbers,:FPL] .> 2.50) .& (households[household_numbers,:FPL] .<= 4)) ; 
				X[:,findall(all_covariates .== "previous_choice_gt400")[1]] = X[:,findall(all_covariates .== "previous_choice")[1]] .* (households[household_numbers,:FPL] .> 4) ; 
				if add_complex_interactions
					X[:,findall(all_covariates .== "previous_choice_0to34")[1]] = X[:,prev_ind] .* (households[household_numbers,:perc_0to17] .+
						households[household_numbers,:perc_18to25] .+ households[household_numbers,:perc_26to34]);
				else
					X[:,findall(all_covariates .== "previous_choice_0to17")[1]] = X[:,prev_ind] .* (households[household_numbers,:perc_0to17]);
					X[:,findall(all_covariates .== "previous_choice_18to34")[1]] = X[:,prev_ind] .* (households[household_numbers,:perc_18to25] .+ households[household_numbers,:perc_26to34]);
				end
				X[:,findall(all_covariates .== "previous_choice_35to54")[1]] = X[:,findall(all_covariates .== "previous_choice")[1]] .* (households[household_numbers,:perc_35to44] .+ households[household_numbers,:perc_45to54]);
				X[:,findall(all_covariates .== "previous_choice_male")[1]] = X[:,findall(all_covariates .== "previous_choice")[1]] .* (households[household_numbers,:perc_male]);
				X[:,findall(all_covariates .== "previous_choice_family")[1]] = X[:,findall(all_covariates .== "previous_choice")[1]] .* (households[household_numbers,:household_size] .> 1);
				X[:,findall(all_covariates .== "previous_choice_Asian")[1]] = X[:,findall(all_covariates .== "previous_choice")[1]] .* (households[household_numbers,:perc_asian]);
				X[:,findall(all_covariates .== "previous_choice_Black")[1]] = X[:,findall(all_covariates .== "previous_choice")[1]] .* (households[household_numbers,:perc_black]);
				X[:,findall(all_covariates .== "previous_choice_Hispanic")[1]] = X[:,findall(all_covariates .== "previous_choice")[1]] .* (households[household_numbers,:perc_hispanic]);
				X[:,findall(all_covariates .== "previous_choice_Other")[1]] = X[:,findall(all_covariates .== "previous_choice")[1]] .* (households[household_numbers,:perc_other]);
			end
			
			if add_complex_interactions
				X[:,findall(all_covariates .== "previous_choice_250to400_0to34")[1]] = X[:,prev_ind] .* ((households[household_numbers,:FPL] .> 2.50) .& (households[household_numbers,:FPL] .<= 4)) .*
					(households[household_numbers,:perc_0to17] .+ households[household_numbers,:perc_18to25] .+ households[household_numbers,:perc_26to34]); 
				X[:,findall(all_covariates .== "previous_choice_gt400_0to34")[1]] = X[:,prev_ind] .* (households[household_numbers,:FPL] .> 4) .*
					(households[household_numbers,:perc_0to17] .+ households[household_numbers,:perc_18to25] .+ households[household_numbers,:perc_26to34]);  
				X[:,findall(all_covariates .== "previous_choice_250to400_35to54")[1]] = X[:,prev_ind] .* ((households[household_numbers,:FPL] .> 2.50) .& (households[household_numbers,:FPL] .<= 4)) .*
					(households[household_numbers,:perc_35to44] .+ households[household_numbers,:perc_45to54]); 
				X[:,findall(all_covariates .== "previous_choice_gt400_35to54")[1]] = X[:,prev_ind] .* (households[household_numbers,:FPL] .> 4) .*
					(households[household_numbers,:perc_35to44] .+ households[household_numbers,:perc_45to54]);  
			end
		end
		
		if add_previous_firm_interactions & !inattention & !add_complex_interactions
			X[:,findall(all_covariates .== "previous_choice_Anthem")[1]] = X[:,prev_ind] .* X[:,findall(all_covariates .== "Anthem")[1]];
			X[:,findall(all_covariates .== "previous_choice_Blue_Shield")[1]] = X[:,prev_ind] .* X[:,findall(all_covariates .== "Blue_Shield")[1]];
			X[:,findall(all_covariates .== "previous_choice_Kaiser")[1]] = X[:,prev_ind] .* X[:,findall(all_covariates .== "Kaiser")[1]];
			X[:,findall(all_covariates .== "previous_choice_Health_Net")[1]] = X[:,prev_ind] .* X[:,findall(all_covariates .== "Health_Net")[1]];			
			X[:,findall(all_covariates .== "previous_choice_HMO")[1]] = X[:,prev_ind] .* X[:,findall(all_covariates .== "HMO")[1]];
			X[:,findall(all_covariates .== "previous_choice_AV")[1]] = X[:,prev_ind] .* X[:,findall(all_covariates .== "AV")[1]];
			if add_silver
				X[:,findall(all_covariates .== "previous_choice_silver")[1]] = X[:,prev_ind] .* X[:,findall(all_covariates .== "Silver")[1]];
			end
		end
	
		# Add intercepts
		if add_intercepts
			if add_complex_interactions
				X[:,findall(all_covariates .== "intercept_year2015")] = (1 .- uninsured_plans) .* (households[household_numbers,:year] .== 2015);
				X[:,findall(all_covariates .== "intercept_year2016")] = (1 .- uninsured_plans) .* (households[household_numbers,:year] .== 2016);
				X[:,findall(all_covariates .== "intercept_year2017")] = (1 .- uninsured_plans) .* (households[household_numbers,:year] .== 2017);
				X[:,findall(all_covariates .== "intercept_year2018")] = (1 .- uninsured_plans) .* (households[household_numbers,:year] .== 2018);
				#X[:,findall(all_covariates .== "intercept_0to34")] = (1 .- uninsured_plans)  .* (households[household_numbers,:perc_0to17] .+
				#	households[household_numbers,:perc_18to25] .+ households[household_numbers,:perc_26to34]);
				#X[:,findall(all_covariates .== "intercept_35to54")] = (1 .- uninsured_plans) .* 
				#	(households[household_numbers,:perc_35to44] .+ households[household_numbers,:perc_45to54]);
				X[:,findall(all_covariates .== "HMO_0to34")[1]] = X[:,findall(all_covariates .== "HMO")[1]] .* (households[household_numbers,:perc_0to17] .+
					households[household_numbers,:perc_18to25] .+ households[household_numbers,:perc_26to34]);
				X[:,findall(all_covariates .== "HMO_35to54")[1]] = X[:,findall(all_covariates .== "HMO")[1]] .* 
					(households[household_numbers,:perc_35to44] .+ households[household_numbers,:perc_45to54]);
			else
				X[:,findall(all_covariates .== "intercept_male")] = (1 .- uninsured_plans) .* (households[household_numbers,:perc_male]);
				X[:,findall(all_covariates .== "intercept_family")] = (1 .- uninsured_plans) .* (households[household_numbers,:household_size] .> 1);
				X[:,findall(all_covariates .== "intercept_0to17")] = (1 .- uninsured_plans)  .* (households[household_numbers,:perc_0to17]);
				X[:,findall(all_covariates .== "intercept_18to34")] = (1 .- uninsured_plans) .* (households[household_numbers,:perc_18to25] .+ households[household_numbers,:perc_26to34]);	
				X[:,findall(all_covariates .== "intercept_35to54")] = (1 .- uninsured_plans) .* (households[household_numbers,:perc_35to44] .+ households[household_numbers,:perc_45to54]);
				X[:,findall(all_covariates .== "intercept_250to400")] = (1 .- uninsured_plans) .* ((households[household_numbers,:FPL] .> 2.5) .& (households[household_numbers,:FPL] .<= 4)) ;						
				X[:,findall(all_covariates .== "intercept_gt400")] = (1 .- uninsured_plans) .* (households[household_numbers,:FPL] .> 4) ; 
				X[:,findall(all_covariates .== "intercept_asian")] = (1 .- uninsured_plans) .* (households[household_numbers,:perc_asian]);
				X[:,findall(all_covariates .== "intercept_black")] = (1 .- uninsured_plans) .* (households[household_numbers,:perc_black]);	
				X[:,findall(all_covariates .== "intercept_hispanic")] = (1 .- uninsured_plans) .* (households[household_numbers,:perc_hispanic]);
				X[:,findall(all_covariates .== "intercept_other")] = (1 .- uninsured_plans) .* (households[household_numbers,:perc_other]);
			end
		end
		
		
		# Add rating areas
		if add_rating_areas
			X[:,findall(all_covariates .== "intercept_ra2")[1]] = (1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 2);
			X[:,findall(all_covariates .== "intercept_ra3")[1]] = (1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 3);
			X[:,findall(all_covariates .== "intercept_ra4")[1]] = (1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 4);
			X[:,findall(all_covariates .== "intercept_ra5")[1]] = (1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 5);
			X[:,findall(all_covariates .== "intercept_ra6")[1]] = (1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 6);
			X[:,findall(all_covariates .== "intercept_ra7")[1]] = (1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 7);
			X[:,findall(all_covariates .== "intercept_ra8")[1]] = (1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 8);
			X[:,findall(all_covariates .== "intercept_ra9")[1]] = (1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 9);
			X[:,findall(all_covariates .== "intercept_ra10")[1]] = (1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 10);
			X[:,findall(all_covariates .== "intercept_ra11")[1]] = (1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 11);
			X[:,findall(all_covariates .== "intercept_ra12")[1]] = (1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 12);
			X[:,findall(all_covariates .== "intercept_ra14")[1]] = (1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 14);
			X[:,findall(all_covariates .== "intercept_ra15")[1]] = (1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 15);
			X[:,findall(all_covariates .== "intercept_ra16")[1]] = (1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 16);
			X[:,findall(all_covariates .== "intercept_ra17")[1]] = (1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 17);
			X[:,findall(all_covariates .== "intercept_ra18")[1]] = (1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 18);
			X[:,findall(all_covariates .== "intercept_ra19")[1]] = (1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 19);
		end	
		
		# Add insurer-market fixed effects
		if add_insurer_market_fe
			X[:,findall(all_covariates .== "intercept_ra2_Anthem")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 2) .*
				X[:,findall(all_covariates .== "Anthem")[1]];
			X[:,findall(all_covariates .== "intercept_ra2_Blue_Shield")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 2) .*
				X[:,findall(all_covariates .== "Blue_Shield")[1]];
			X[:,findall(all_covariates .== "intercept_ra2_Health_Net")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 2) .*
				X[:,findall(all_covariates .== "Health_Net")[1]];	
			X[:,findall(all_covariates .== "intercept_ra2_Kaiser")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 2) .*
				X[:,findall(all_covariates .== "Kaiser")[1]];

			X[:,findall(all_covariates .== "intercept_ra3_Anthem")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 3) .*
				X[:,findall(all_covariates .== "Anthem")[1]];
			X[:,findall(all_covariates .== "intercept_ra3_Blue_Shield")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 3) .*
				X[:,findall(all_covariates .== "Blue_Shield")[1]];
			X[:,findall(all_covariates .== "intercept_ra3_Kaiser")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 3) .*
				X[:,findall(all_covariates .== "Kaiser")[1]];
			
			X[:,findall(all_covariates .== "intercept_ra4_Anthem")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 4) .*
				X[:,findall(all_covariates .== "Anthem")[1]];
			X[:,findall(all_covariates .== "intercept_ra4_Blue_Shield")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 4) .*
				X[:,findall(all_covariates .== "Blue_Shield")[1]];
			X[:,findall(all_covariates .== "intercept_ra4_Health_Net")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 4) .*
				X[:,findall(all_covariates .== "Health_Net")[1]];	
			X[:,findall(all_covariates .== "intercept_ra4_Kaiser")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 4) .*
				X[:,findall(all_covariates .== "Kaiser")[1]];

			X[:,findall(all_covariates .== "intercept_ra5_Anthem")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 5) .*
				X[:,findall(all_covariates .== "Anthem")[1]];
			X[:,findall(all_covariates .== "intercept_ra5_Blue_Shield")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 5) .*
				X[:,findall(all_covariates .== "Blue_Shield")[1]];
			X[:,findall(all_covariates .== "intercept_ra5_Health_Net")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 5) .*
				X[:,findall(all_covariates .== "Health_Net")[1]];	
					
			X[:,findall(all_covariates .== "intercept_ra6_Anthem")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 6) .*
				X[:,findall(all_covariates .== "Anthem")[1]];
			X[:,findall(all_covariates .== "intercept_ra6_Blue_Shield")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 6) .*
				X[:,findall(all_covariates .== "Blue_Shield")[1]];
			
			X[:,findall(all_covariates .== "intercept_ra7_Anthem")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 7) .*
				X[:,findall(all_covariates .== "Anthem")[1]];
			X[:,findall(all_covariates .== "intercept_ra7_Blue_Shield")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 7) .*
				X[:,findall(all_covariates .== "Blue_Shield")[1]];
			X[:,findall(all_covariates .== "intercept_ra7_Health_Net")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 7) .*
				X[:,findall(all_covariates .== "Health_Net")[1]];	
			X[:,findall(all_covariates .== "intercept_ra7_Kaiser")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 7) .*
				X[:,findall(all_covariates .== "Kaiser")[1]];	
						
			X[:,findall(all_covariates .== "intercept_ra8_Anthem")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 8) .*
				X[:,findall(all_covariates .== "Anthem")[1]];
			X[:,findall(all_covariates .== "intercept_ra8_Blue_Shield")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 8) .*
				X[:,findall(all_covariates .== "Blue_Shield")[1]];
			X[:,findall(all_covariates .== "intercept_ra8_Health_Net")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 8) .*
				X[:,findall(all_covariates .== "Health_Net")[1]];	
			X[:,findall(all_covariates .== "intercept_ra8_Kaiser")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 8) .*
				X[:,findall(all_covariates .== "Kaiser")[1]];
			
			X[:,findall(all_covariates .== "intercept_ra9_Anthem")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 9) .*
				X[:,findall(all_covariates .== "Anthem")[1]];
			X[:,findall(all_covariates .== "intercept_ra9_Blue_Shield")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 9) .*
				X[:,findall(all_covariates .== "Blue_Shield")[1]];
			
			X[:,findall(all_covariates .== "intercept_ra10_Anthem")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 10) .*
				X[:,findall(all_covariates .== "Anthem")[1]];
			X[:,findall(all_covariates .== "intercept_ra10_Blue_Shield")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 10) .*
				X[:,findall(all_covariates .== "Blue_Shield")[1]];
			X[:,findall(all_covariates .== "intercept_ra10_Health_Net")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 10) .*
				X[:,findall(all_covariates .== "Health_Net")[1]];		
			
			X[:,findall(all_covariates .== "intercept_ra11_Anthem")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 11) .*
				X[:,findall(all_covariates .== "Anthem")[1]];
			X[:,findall(all_covariates .== "intercept_ra11_Blue_Shield")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 11) .*
				X[:,findall(all_covariates .== "Blue_Shield")[1]];
			
			X[:,findall(all_covariates .== "intercept_ra12_Anthem")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 12) .*
				X[:,findall(all_covariates .== "Anthem")[1]];
			X[:,findall(all_covariates .== "intercept_ra12_Blue_Shield")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 12) .*
				X[:,findall(all_covariates .== "Blue_Shield")[1]];
			
			X[:,findall(all_covariates .== "intercept_ra14_Anthem")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 14) .*
				X[:,findall(all_covariates .== "Anthem")[1]];
			X[:,findall(all_covariates .== "intercept_ra14_Blue_Shield")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 14) .*
				X[:,findall(all_covariates .== "Blue_Shield")[1]];
			X[:,findall(all_covariates .== "intercept_ra14_Health_Net")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 14) .*
				X[:,findall(all_covariates .== "Health_Net")[1]];
			
			X[:,findall(all_covariates .== "intercept_ra15_Anthem")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 15) .*
				X[:,findall(all_covariates .== "Anthem")[1]];
			X[:,findall(all_covariates .== "intercept_ra15_Blue_Shield")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 15) .*
				X[:,findall(all_covariates .== "Blue_Shield")[1]];
			X[:,findall(all_covariates .== "intercept_ra15_Kaiser")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 15) .*
				X[:,findall(all_covariates .== "Kaiser")[1]];
			X[:,findall(all_covariates .== "intercept_ra15_Health_Net")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 15) .*
				X[:,findall(all_covariates .== "Health_Net")[1]];
			
			X[:,findall(all_covariates .== "intercept_ra16_Anthem")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 16) .*
				X[:,findall(all_covariates .== "Anthem")[1]];
			X[:,findall(all_covariates .== "intercept_ra16_Blue_Shield")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 16) .*
				X[:,findall(all_covariates .== "Blue_Shield")[1]];							
			X[:,findall(all_covariates .== "intercept_ra16_Kaiser")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 16) .*
				X[:,findall(all_covariates .== "Kaiser")[1]];
			X[:,findall(all_covariates .== "intercept_ra16_Health_Net")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 16) .*
				X[:,findall(all_covariates .== "Health_Net")[1]];
				
			X[:,findall(all_covariates .== "intercept_ra17_Anthem")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 17) .*
				X[:,findall(all_covariates .== "Anthem")[1]];
			X[:,findall(all_covariates .== "intercept_ra17_Blue_Shield")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 17) .*
				X[:,findall(all_covariates .== "Blue_Shield")[1]];
			X[:,findall(all_covariates .== "intercept_ra17_Health_Net")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 17) .*
				X[:,findall(all_covariates .== "Health_Net")[1]];
			X[:,findall(all_covariates .== "intercept_ra17_Kaiser")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 17) .*
				X[:,findall(all_covariates .== "Kaiser")[1]];
			
			X[:,findall(all_covariates .== "intercept_ra18_Anthem")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 18) .*
				X[:,findall(all_covariates .== "Anthem")[1]];
			X[:,findall(all_covariates .== "intercept_ra18_Blue_Shield")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 18) .*
				X[:,findall(all_covariates .== "Blue_Shield")[1]];
			X[:,findall(all_covariates .== "intercept_ra18_Health_Net")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 18) .*
				X[:,findall(all_covariates .== "Health_Net")[1]];			
			X[:,findall(all_covariates .== "intercept_ra18_Kaiser")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 18) .*
				X[:,findall(all_covariates .== "Kaiser")[1]];
				
			#X[:,findall(all_covariates .== "intercept_ra19_Anthem")[1]] =
			#	(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 19) .*
			#	X[:,findall(all_covariates .== "Anthem")[1]];
			X[:,findall(all_covariates .== "intercept_ra19_Blue_Shield")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 19) .*
				X[:,findall(all_covariates .== "Blue_Shield")[1]];
			X[:,findall(all_covariates .== "intercept_ra19_Kaiser")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 19) .*
				X[:,findall(all_covariates .== "Kaiser")[1]];
			X[:,findall(all_covariates .== "intercept_ra19_Health_Net")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 19) .*
				X[:,findall(all_covariates .== "Health_Net")[1]];
		end
		
		# Residual from reduced form stage
		if control_function
			X[:,findall(all_covariates .== "residual")[1]] = convert(Array,data[!,:residual]);
		end
		
		# Random Parameter Variables
		if random_parameters
			X[:,findall(all_covariates .== "intercept_random")[1]] = X[:,findall(all_covariates .== "intercept")[1]];
			X[:,findall(all_covariates .== "Anthem_random")[1]] = X[:,findall(all_covariates .== "Anthem")[1]];
			X[:,findall(all_covariates .== "Blue_Shield_random")[1]] = X[:,findall(all_covariates .== "Blue_Shield")[1]];
			X[:,findall(all_covariates .== "Kaiser_random")[1]] = X[:,findall(all_covariates .== "Kaiser")[1]];
			X[:,findall(all_covariates .== "Health_Net_random")[1]] = X[:,findall(all_covariates .== "Health_Net")[1]];
			if !add_av_interactions
				X[:,findall(all_covariates .== "AV_random")[1]] = X[:,findall(all_covariates .== "AV")[1]];
			end
			if add_silver
				X[:,findall(all_covariates .== "Silver_random")[1]] = X[:,findall(all_covariates .== "Silver")[1]];
			end
		end
		
		return(X)
	end	
	
##### Output objects

	fields = ["Total_Enrollment";"Total_Renewals";"MSE";"AV";"AVnocsr";"Bronze";"Silver";"Gold";"Platinum";
				"Anthem";"Blue_Shield";"Health_Net";"Kaiser";"Other_Insurer";
				"Switching Rate";"Avg Premium"];
	cols = ["Observed","Simulated"];
	output = zeros(length(fields),length(cols));


	#### Term for control function approach	
	if control_function
		eta = Array(randn(length(choices),ns));		
	end	
		
	#### Terms for Random Parameters
	if random_parameters 
		
		# Assign a household reference number
		household_ref_numbers = zeros(maximum(households[!,:household_id]));
		household_ref_numbers[unique(households[!,:household_id])] = collect(1:1:length(unique(households[!,:household_id])));
		households[!,:ref_number] = Int64.(household_ref_numbers[households[!,:household_id]]);
		
		V = zeros(length(household_numbers),length(random_variables),ns);
		for k in 1:length(random_variables)
			Random.seed!(random_parameter_product_indices[k]);
			Vi = randn(length(unique(households[!,:ref_number])),ns);
			Vi = Vi[households[!,:ref_number],:];
			V[:,k,:] = Vi[household_numbers,:];
		end
	else
		V = 0;
	end
	
	##### Create Covariate Matrix
	X = form_X_matrix(all_covariates,data,households,choices,previous_choices,weights,uninsured_plans,household_numbers,V);
	data[!,:unique_id] = string.(data[!,:household_number],"_",data[!,:plan_name]);
	probabilities = zeros(size(data)[1]);
	previous_probabilities = zeros(size(data)[1]);
	data[!,:prob] = probabilities;
	
	##### Create previous choice links
	data[!,:previous_choice_index] .= 0;
	data[!,:household_id] = households[household_numbers,:household_id];
	data[!,:unique_id_year] = string.(data[!,:household_id],"_",data[!,:plan_name],"_",data[!,:year]);
	for yr in 2015:2018
		year_data_indices = findall(data[!,:year] .== yr);
		prev_year_data_indices = indexin(string.(data[year_data_indices,:household_id],"_",data[year_data_indices,:plan_name],"_",data[year_data_indices,:year] .- 1),data[!,:unique_id_year]);
		prev_year_data_indices[prev_year_data_indices .== nothing] .= 0;
		data[year_data_indices,:previous_choice_index] = Int64.(prev_year_data_indices); 
	end		
	
	##### Output Processing
		

		Random.seed!(1);
		#households[households[!,:plan_name] .== "NA",:plan_name] .= "Uninsured";
		households[!,:simulated_plan_choice] = copy(convert(Array,households[!,:plan_name]));
		households[!,:simulated_plan_choice] .= "TBD";
		households[!,:simulated_previous_plan_choice] .= "TBD";
		households[!,:unique_id] = string.(households[!,:household_id],"_",households[!,:year]);
		
		for yr in 2014:2018
		
			Xupdated = copy(X);
		
			# Get indices
			households_indices = findall(households[!,:year] .== yr);
			household_data_indices = findall(in(households_indices),household_numbers);
		
			# Update previous choice variables for previously enrolled households (in data and not uninsured)
			if (yr > 2014) & forward_flag
				#previous_households_indices = intersect(findall(households[!,:year] .== yr - 1),findall(households[!,:simulated_plan_choice] .!= "Uninsured"));
				
				
				previous_households_indices = findall(households[!,:year] .== yr - 1);
				previous_households_unique_ids = households[previous_households_indices,:unique_id];
				previous_households_plans = households[previous_households_indices,:simulated_plan_choice];
				
				# Determine renewing households and new households (could be uninsured or not in market)
				renewing_households_indices = households_indices[findall(in(previous_households_unique_ids),string.(households[households_indices,:household_id],"_",yr-1))];
				new_households_indices = setdiff(households_indices,renewing_households_indices);
			
				households[renewing_households_indices,:simulated_previous_plan_choice] = 
					convert(Array,households[indexin(string.(households[renewing_households_indices,:household_id],"_",yr-1),households[!,:unique_id]),:simulated_plan_choice]);
				households[new_households_indices,:simulated_previous_plan_choice] .= "NA";
			
				# Update previous choice variable for renewing households
				prev_ind = findall(all_covariates .== "previous_choice")[1];
				previous_choices_data_ids = string.(renewing_households_indices,"_",households[renewing_households_indices,:simulated_previous_plan_choice]);
				Xupdated[household_data_indices,prev_ind] .= 0;	
				Xupdated[findall(in(previous_choices_data_ids),data[!,:unique_id]),prev_ind] .= 1;
				Xupdated[uninsured_plans .== 1,prev_ind] .== 0;	

				# Update Interaction Terms
				Xupdated[household_data_indices,findall(all_covariates .== "previous_choice_250to400")[1]] = Xupdated[household_data_indices,prev_ind] .* ((households[household_numbers[household_data_indices],:FPL] .> 2.50) .& (households[household_numbers[household_data_indices],:FPL] .<= 4)) ; 
				Xupdated[household_data_indices,findall(all_covariates .== "previous_choice_gt400")[1]] = Xupdated[household_data_indices,prev_ind] .* (households[household_numbers[household_data_indices],:FPL] .> 4) ; 
				if add_complex_interactions
					Xupdated[household_data_indices,findall(all_covariates .== "previous_choice_0to34")[1]] = Xupdated[household_data_indices,prev_ind] .* (households[household_numbers[household_data_indices],:perc_0to17] .+
						households[household_numbers[household_data_indices],:perc_18to25] .+ households[household_numbers[household_data_indices],:perc_26to34]);
				else
					Xupdated[household_data_indices,findall(all_covariates .== "previous_choice_0to17")[1]] = Xupdated[household_data_indices,prev_ind] .* (households[household_numbers[household_data_indices],:perc_0to17]);
					Xupdated[household_data_indices,findall(all_covariates .== "previous_choice_18to34")[1]] = Xupdated[household_data_indices,prev_ind] .* (households[household_numbers[household_data_indices],:perc_18to25] .+ households[household_numbers[household_data_indices],:perc_26to34]);
				end
				Xupdated[household_data_indices,findall(all_covariates .== "previous_choice_35to54")[1]] = Xupdated[household_data_indices,prev_ind] .* (households[household_numbers[household_data_indices],:perc_35to44] .+ households[household_numbers[household_data_indices],:perc_45to54]);
				Xupdated[household_data_indices,findall(all_covariates .== "previous_choice_male")[1]] = Xupdated[household_data_indices,prev_ind] .* (households[household_numbers[household_data_indices],:perc_male]);
				Xupdated[household_data_indices,findall(all_covariates .== "previous_choice_family")[1]] = Xupdated[household_data_indices,prev_ind] .* (households[household_numbers[household_data_indices],:household_size] .> 1);
				Xupdated[household_data_indices,findall(all_covariates .== "previous_choice_Asian")[1]] = Xupdated[household_data_indices,prev_ind] .* (households[household_numbers[household_data_indices],:perc_asian]);
				Xupdated[household_data_indices,findall(all_covariates .== "previous_choice_Black")[1]] = Xupdated[household_data_indices,prev_ind] .* (households[household_numbers[household_data_indices],:perc_black]);
				Xupdated[household_data_indices,findall(all_covariates .== "previous_choice_Hispanic")[1]] = Xupdated[household_data_indices,prev_ind] .* (households[household_numbers[household_data_indices],:perc_hispanic]);
				Xupdated[household_data_indices,findall(all_covariates .== "previous_choice_Other")[1]] = Xupdated[household_data_indices,prev_ind] .* (households[household_numbers[household_data_indices],:perc_other]);
				
				
				if add_complex_interactions
					Xupdated[household_data_indices,findall(all_covariates .== "previous_choice_250to400_0to34")[1]] = Xupdated[household_data_indices,prev_ind] .* ((households[household_numbers[household_data_indices],:FPL] .> 2.50) .& (households[household_numbers[household_data_indices],:FPL] .<= 4)) .*
						(households[household_numbers[household_data_indices],:perc_0to17] .+ households[household_numbers[household_data_indices],:perc_18to25] .+ households[household_numbers[household_data_indices],:perc_26to34]); 
					Xupdated[household_data_indices,findall(all_covariates .== "previous_choice_gt400_0to34")[1]] = Xupdated[household_data_indices,prev_ind] .* (households[household_numbers[household_data_indices],:FPL] .> 4) .*
						(households[household_numbers[household_data_indices],:perc_0to17] .+ households[household_numbers[household_data_indices],:perc_18to25] .+ households[household_numbers[household_data_indices],:perc_26to34]);  
					Xupdated[household_data_indices,findall(all_covariates .== "previous_choice_250to400_35to54")[1]] = Xupdated[household_data_indices,prev_ind] .* ((households[household_numbers[household_data_indices],:FPL] .> 2.50) .& (households[household_numbers[household_data_indices],:FPL] .<= 4)) .*
						(households[household_numbers[household_data_indices],:perc_35to44] .+ households[household_numbers[household_data_indices],:perc_45to54]); 
					Xupdated[household_data_indices,findall(all_covariates .== "previous_choice_gt400_35to54")[1]] = Xupdated[household_data_indices,prev_ind] .* (households[household_numbers[household_data_indices],:FPL] .> 4) .*
						(households[household_numbers[household_data_indices],:perc_35to44] .+ households[household_numbers[household_data_indices],:perc_45to54]);  
				end
			end
			
			# Simulate choice probabilities
			probabilities[household_data_indices] = get_choice_probabilities(min_beta,Xupdated,V,choices,household_numbers,uninsured_plans)[household_data_indices];
			
			
			# Now use the choice probabilities to simulate a choice
			data[household_data_indices,:prob] = probabilities[household_data_indices];
			simulated_plan_choices = combine(groupby(data[household_data_indices,:], :household_number), [:plan_name, :prob] => ((x, y) -> sample(x, Weights(y))) => :simulated_plan_choice);
			households[households_indices,:simulated_plan_choice] = simulated_plan_choices[!,:simulated_plan_choice];

		end
		
		# Switching Stats
		data[!,:previous_prob] .= zeros(size(data)[1]);
		data[findall(data[!,:previous_choice_index] .> 0),:previous_prob] = data[data[findall(data[!,:previous_choice_index] .> 0),:previous_choice_index],:prob];
		data[!,:prob_diff] = abs.(data[!,:prob] .- data[!,:previous_prob]);
		
		household_previous_prob = combine(groupby(data,:household_number), :previous_prob => sum)[!,:previous_prob_sum];
		household_prob = combine(groupby(data,:household_number), :prob => sum)[!,:prob_sum];
		
		household_previous_prob_nouni = combine(groupby(data[uninsured_plans .== 0,:],:household_number), :previous_prob => sum)[!,:previous_prob_sum];
		household_prob_nouni = combine(groupby(data[uninsured_plans .== 0,:],:household_number), :prob => sum)[!,:prob_sum];
		
		enrolled_prev_year = findall(household_previous_prob .> 0);
		enrolled_prev_year_data_indices = intersect(findall(in(enrolled_prev_year),household_numbers),findall(uninsured_plans .== 0));
		
		prob_not_offered = household_prob .- household_previous_prob; # note this will be 1 for new enrollees, but this doesn't matter to calculations below
		prob_uni_switch = abs.(household_prob_nouni .- household_previous_prob_nouni); # this could be negative of probability of uninsured goes up
		
		household_switching_rates = max.((combine(groupby(data[enrolled_prev_year_data_indices,:],:household_number), :prob_diff => sum)[!,:prob_diff_sum] .- prob_not_offered[enrolled_prev_year] .- prob_uni_switch[enrolled_prev_year])./2,0);
		
		
		
		# Extract simulated choice and previous choices
		simulated_choice_data_ids = string.(collect(1:1:size(households)[1]),"_",households[!,:simulated_plan_choice]);
		simulated_choices = zeros(length(choices));
		simulated_choices[findall(in(simulated_choice_data_ids),data[!,:unique_id])] .= 1;
		
		simulated_previous_choice_data_ids = string.(collect(1:1:size(households)[1]),"_",households[!,:simulated_previous_plan_choice]);
		simulated_previous_choices = zeros(length(choices));
		simulated_previous_choices[findall(in(simulated_previous_choice_data_ids),data[!,:unique_id])] .= 1;
		simulated_previous_choices[uninsured_plans .== 1] .= 0;
		
		# Get Switching indices (observed and simulated)
		previous_indices = findall(previous_choices .== 1); # people who renewed coverage
		previous_households = household_numbers[previous_indices]; 
		previously_enrolled_indices = findall(in(previous_households),household_numbers);
		renewal_indices = intersect(findall(uninsured_plans .== 0),previously_enrolled_indices); # simulated people who renewed coverage
		same_indices = intersect(renewal_indices,previous_indices);
		
		simulated_previous_indices = findall(simulated_previous_choices .== 1); # people who renewed coverage
		simulated_previous_households = household_numbers[simulated_previous_indices]; 
		simulated_previously_enrolled_indices = findall(in(simulated_previous_households),household_numbers);
		simulated_renewal_indices = intersect(findall(uninsured_plans .== 0),simulated_previously_enrolled_indices); # simulated people who renewed coverage
		simulated_same_indices = intersect(simulated_renewal_indices,simulated_previous_indices);
		
		choice_weights = choices .* weights;
		simulated_choice_weights = simulated_choices .* weights;			
		println("Done")
		
		avnocsr = X[:,findall(all_covariates .== "AV")[1]];
		avnocsr[avnocsr .== 0.73] .= 0.7;
		avnocsr[avnocsr .== 0.87] .= 0.7;
		avnocsr[avnocsr .== 0.94] .= 0.7;
		
		# Actual Output
		output[findall(fields .== "Total_Enrollment")[1],findall(cols .== "Observed")[1]] += sum(choice_weights[uninsured_plans .== 0]);
		output[findall(fields .== "Total_Renewals")[1],findall(cols .== "Observed")[1]] += sum(choice_weights[renewal_indices]);

		output[findall(fields .== "AV")[1],findall(cols .== "Observed")[1]] = 
			sum(choice_weights[uninsured_plans .== 0] .* X[uninsured_plans .== 0,findall(all_covariates .== "AV")[1]])/
				sum(choice_weights[uninsured_plans .== 0]);
		output[findall(fields .== "AVnocsr")[1],findall(cols .== "Observed")[1]] = 
			sum(choice_weights[uninsured_plans .== 0] .* avnocsr[uninsured_plans .== 0])/sum(choice_weights[uninsured_plans .== 0]);
		output[findall(fields .== "Bronze")[1],findall(cols .== "Observed")[1]] += sum(choice_weights[findall(data[!,:av] .== 0.6)]);
		output[findall(fields .== "Gold")[1],findall(cols .== "Observed")[1]] += sum(choice_weights[findall(data[!,:av] .== 0.8)]);
		output[findall(fields .== "Platinum")[1],findall(cols .== "Observed")[1]] += sum(choice_weights[findall(data[!,:av] .== 0.9)]);
		output[findall(fields .== "Silver")[1],findall(cols .== "Observed")[1]] += 
			sum(choice_weights[cat(dims=1,findall(data[!,:av] .== 0.7),findall(data[!,:av] .== 0.73),
				findall(data[!,:av] .== 0.87),findall(data[!,:av] .== 0.94))]);

		output[findall(fields .== "Anthem")[1],findall(cols .== "Observed")[1]] += 
			sum(choice_weights .* X[:,findall(all_covariates .== "Anthem")[1]]);
		output[findall(fields .== "Blue_Shield")[1],findall(cols .== "Observed")[1]] += 
			sum(choice_weights .* X[:,findall(all_covariates .== "Blue_Shield")[1]]);
		output[findall(fields .== "Health_Net")[1],findall(cols .== "Observed")[1]] += 
			sum(choice_weights .* X[:,findall(all_covariates .== "Health_Net")[1]]);
		output[findall(fields .== "Kaiser")[1],findall(cols .== "Observed")[1]] += 
			sum(choice_weights .* X[:,findall(all_covariates .== "Kaiser")[1]]);
		output[findall(fields .== "Other_Insurer")[1],findall(cols .== "Observed")[1]] += 
			sum(choice_weights[uninsured_plans .== 0]) - sum(choice_weights .* 
				(X[:,findall(all_covariates .== "Anthem")[1]] .+ X[:,findall(all_covariates .== "Blue_Shield")[1]] .+ 
				X[:,findall(all_covariates .== "Health_Net")[1]] .+ X[:,findall(all_covariates .== "Kaiser")[1]]));
		
		prem_ind = findall(all_covariates .== "Premium")[1];
		output[findall(fields .== "Avg Premium")[1],findall(cols .== "Observed")[1]] += sum(choice_weights[uninsured_plans .== 0] .* 
			(X[uninsured_plans .== 0,prem_ind] .+ data[uninsured_plans .== 0,:penalty] ./ weights[uninsured_plans .== 0]));

		output[findall(fields .== "Switching Rate")[1],findall(cols .== "Observed")[1]] += sum(choice_weights[same_indices]);
		#output[findall(fields .== "Switching Rate")[1],findall(cols .== "Observed")[1]] += 
		#	sum(choice_weights)/sum(households[!,:weight]) * (sum(choice_weights) - sum(choice_weights[same_indices]))/(2 * sum(choice_weights));
			
			
			
		# Simulated Output
		output[findall(fields .== "Total_Enrollment")[1],findall(cols .== "Simulated")[1]] += sum(simulated_choice_weights[uninsured_plans .== 0]);
		output[findall(fields .== "Total_Renewals")[1],findall(cols .== "Simulated")[1]] += sum(simulated_choice_weights[renewal_indices]);
		output[findall(fields .== "MSE")[1],findall(cols .== "Simulated")[1]] = 
			#1/sum(weights[choices .== 1]) .* sum(((choices[choices .== 1] .- simulated_choices[choices .== 1]).^2) .* weights[choices .== 1]);
			1/sum(weights) .* sum(((choices .- simulated_choices).^2) .* weights);
	
		output[findall(fields .== "AV")[1],findall(cols .== "Simulated")[1]] = 
			sum(simulated_choice_weights[uninsured_plans .== 0] .* X[uninsured_plans .== 0,findall(all_covariates .== "AV")[1]])/
				sum(simulated_choice_weights[uninsured_plans .== 0]);
		output[findall(fields .== "AVnocsr")[1],findall(cols .== "Simulated")[1]] = 
			sum(simulated_choice_weights[uninsured_plans .== 0] .* avnocsr[uninsured_plans .== 0])/sum(simulated_choice_weights[uninsured_plans .== 0]);
		output[findall(fields .== "Bronze")[1],findall(cols .== "Simulated")[1]] += sum(simulated_choice_weights[findall(data[!,:av] .== 0.6)]);
		output[findall(fields .== "Gold")[1],findall(cols .== "Simulated")[1]] += sum(simulated_choice_weights[findall(data[!,:av] .== 0.8)]);
		output[findall(fields .== "Platinum")[1],findall(cols .== "Simulated")[1]] += sum(simulated_choice_weights[findall(data[!,:av] .== 0.9)]);
		output[findall(fields .== "Silver")[1],findall(cols .== "Simulated")[1]] += 
			sum(simulated_choice_weights[cat(dims=1,findall(data[!,:av] .== 0.7),findall(data[!,:av] .== 0.73),
				findall(data[!,:av] .== 0.87),findall(data[!,:av] .== 0.94))]);
		output[findall(in(["Bronze","Silver","Gold","Platinum"]),fields),findall(cols .== "Simulated")[1]] += get_dominated_enrollment(simulated_choice_weights,X,data,households,household_numbers);
		
		output[findall(fields .== "Anthem")[1],findall(cols .== "Simulated")[1]] += 
			sum(simulated_choice_weights .* X[:,findall(all_covariates .== "Anthem")[1]]);
		output[findall(fields .== "Blue_Shield")[1],findall(cols .== "Simulated")[1]] += 
			sum(simulated_choice_weights .* X[:,findall(all_covariates .== "Blue_Shield")[1]]);
		output[findall(fields .== "Health_Net")[1],findall(cols .== "Simulated")[1]] += 
			sum(simulated_choice_weights .* X[:,findall(all_covariates .== "Health_Net")[1]]);
		output[findall(fields .== "Kaiser")[1],findall(cols .== "Simulated")[1]] += 
			sum(simulated_choice_weights .* X[:,findall(all_covariates .== "Kaiser")[1]]);
		output[findall(fields .== "Other_Insurer")[1],findall(cols .== "Simulated")[1]] += 
			sum(simulated_choice_weights[uninsured_plans .== 0]) - sum(simulated_choice_weights .* 
				(X[:,findall(all_covariates .== "Anthem")[1]] .+ X[:,findall(all_covariates .== "Blue_Shield")[1]] .+ 
				X[:,findall(all_covariates .== "Health_Net")[1]] .+ X[:,findall(all_covariates .== "Kaiser")[1]]));
				
		output[findall(fields .== "Avg Premium")[1],findall(cols .== "Simulated")[1]] += sum(simulated_choice_weights[uninsured_plans .== 0] .* 
			(X[uninsured_plans .== 0,prem_ind] .+ data[uninsured_plans .== 0,:penalty] ./ weights[uninsured_plans .== 0]));

		output[findall(fields .== "Switching Rate")[1],findall(cols .== "Simulated")[1]] += sum(simulated_choice_weights[simulated_same_indices]);
		#output[findall(fields .== "Switching Rate")[1],findall(cols .== "Simulated")[1]] += 
		#	sum(simulated_choice_weights)/sum(households[!,:weight]) * (sum(simulated_choice_weights) - sum(simulated_choice_weights[same_indices]))/(2 * sum(simulated_choice_weights));
	
# Compute Shares

	fit_output = DataFrame();
	fit_output[!,:Measure] = fields;
	fit_output[!,:Observed] = output[:,findall(cols .== "Observed")[1]];
	fit_output[!,:Simulated] = output[:,findall(cols .== "Simulated")[1]];

	share_fields = findall(in(["Bronze","Silver","Gold","Platinum",
		"Anthem","Blue_Shield","Health_Net","Kaiser","Other_Insurer","Avg Premium"]),fields);
	fit_output[share_fields,:Observed] = output[share_fields,findall(cols .== "Observed")[1]] ./ 
		output[findall(fields .== "Total_Enrollment")[1],findall(cols .== "Observed")[1]];
	fit_output[share_fields,:Simulated] = output[share_fields,findall(cols .== "Simulated")[1]] ./ 
		output[findall(fields .== "Total_Enrollment")[1],findall(cols .== "Simulated")[1]];

	switching_rate_field = findall(fields .== "Switching Rate")[1];
	fit_output[switching_rate_field,:Observed] =  
		(output[findall(fields .== "Total_Renewals")[1],findall(cols .== "Observed")[1]] .- 
			output[switching_rate_field,findall(cols .== "Observed")[1]]) ./ 
		output[findall(fields .== "Total_Renewals")[1],findall(cols .== "Observed")[1]];
		
	fit_output[switching_rate_field,:Simulated] = 
		(output[findall(fields .== "Total_Renewals")[1],findall(cols .== "Simulated")[1]] .- 
			output[switching_rate_field,findall(cols .== "Simulated")[1]]) ./ 
		output[findall(fields .== "Total_Renewals")[1],findall(cols .== "Simulated")[1]];
	fit_output[switching_rate_field,:Simulated] = sum(household_switching_rates .* households[enrolled_prev_year,:weight])/sum(households[enrolled_prev_year,:weight]);	
	#fit_output[switching_rate_field,:Observed] = output[switching_rate_field,findall(cols .== "Observed")[1]];
	#fit_output[switching_rate_field,:Simulated] = output[switching_rate_field,findall(cols .== "Simulated")[1]];
	
	
	if add_av_interactions & add_insurer_market_fe & sequence_flag & random_parameters & control_function
		CSV.write("model_fit_output_imfe_seq_cf_random_av.csv",fit_output);
	elseif add_av_interactions & sequence_flag & random_parameters & control_function
		CSV.write("model_fit_output_seq_cf_random_av.csv",fit_output);
	elseif add_av_interactions & add_insurer_market_fe & random_parameters & control_function
		CSV.write("model_fit_output_imfe_cf_random_av.csv",fit_output);
	elseif add_insurer_market_fe & sequence_flag & random_parameters & control_function
		CSV.write("model_fit_output_imfe_seq_cf_random.csv",fit_output);
	elseif add_insurer_market_fe & random_parameters & control_function
		CSV.write("model_fit_output_imfe_cf_random.csv",fit_output);
	elseif add_av_interactions & sequence_flag & random_parameters & control_function
		CSV.write("model_fit_output_av_seq_cf_random.csv",fit_output);
	elseif add_av_interactions & random_parameters & control_function
		CSV.write("model_fit_output_av_cf_random.csv",fit_output);
	elseif inattention
		CSV.write("model_fit_output_inattention.csv",fit_output);
	elseif sequence_flag & random_parameters & control_function & add_silver
		CSV.write("model_fit_output_seq_cf_random_silver.csv",fit_output);
	elseif sequence_flag & random_parameters & control_function
		CSV.write("model_fit_output_seq_cf_random.csv",fit_output);
	elseif random_parameters & sequence_flag
		CSV.write("model_fit_output_seq_random.csv",fit_output);
	elseif random_parameters & control_function & add_silver
		CSV.write("model_fit_output_cf_random_silver.csv",fit_output);
	elseif random_parameters & control_function
		CSV.write("model_fit_output_cf_random.csv",fit_output);
	elseif random_parameters
		CSV.write("model_fit_output_random.csv",fit_output);
	elseif nested
		CSV.write("model_fit_output_nested.csv",fit_output);	
	else
		CSV.write("model_fit_output_baselogit.csv",fit_output);	
	end	
	